CRStest <- function(stats, alpha = .05, random = FALSE) {
	q <- length(stats)
	G <- t(unname(expand.grid(rep(list(c(1, -1)), q))))
	ts <- apply(G * stats, 2, function(x) mean(x)/sd(x))
	M <- length(ts)
	k <- ceiling(M * (1 - alpha))
	t.orig <- ts[1]
	ts <- sort(ts)
	as.numeric(t.orig > ts[k])
}


args <- commandArgs(TRUE)
j <- eval( parse(text=args[1]) )
output <- paste("did-crs-", j, ".txt", sep="")

set.seed(48109)

q1 <- 6
q0 <- 6
t0 <- 10
t1 <- 10
rho <- .5
gamma <- .8
beta <- c(1,1,1)
theta <- 1
n.deltas <- 4
deltas <- seq(0, 3, length.out = n.deltas)

sims <- 10000
q <- q0 + q1
n <- t0 + t1
brn <- 500
nbrn <- brn + n 


sigv  <- c(rep(1, q - j), rep(20, j))
sigw  <- c(rep(1, q - j), rep(20, j))
sigx1 <- c(rep(1, q - j), rep(20, j))
sigx2 <- c(rep(1, q - j), rep(20, j))
sigx3 <- c(rep(1, q - j), rep(20, j))

id <- c(rep(1:q1, rep(n, q1)), rep(1:q1, rep(n, q1)))
treat <- rep(0, q*n)
treat[1:(q1*n)] <- 1
iv <- rep(c(rep(0, t0), rep(1, t1)), q)
ivtreat <- iv*treat

totres <- rep(NA, n.deltas)

for(h in 1:n.deltas) {
	res <- rep(NA, sims)
	delta <- deltas[h]

	for(i in 1:sims) {
		u <- matrix(NA, nbrn, q)
		u[1, ] <- 0
		for(t in 1:(nbrn-1)) {
			v <- rnorm(q, 0, sigv)
			u[t+1, ] <- rho*u[t, ] + v
		}
		u <- u[(brn+1):nbrn, ]
		u <- c(u)
	
		x1 <- matrix(NA, n, q)
		x2 <- matrix(NA, n, q)
		x3 <- matrix(NA, n, q)
		mat.ivtreat <- matrix(ivtreat, n, q)
		for(k in 1:q) { 
			x1[, k] <- gamma*mat.ivtreat[, k] + rnorm(n, 0, sigx1[k])
			x2[, k] <- rnorm(n, 0, sigx2[k])
			x3[, k] <- rnorm(n, 0, sigx3[k])
			}
		x1 <- c(x1)
		x2 <- c(x2)
		x3 <- c(x3)
		
		y <- 1 + theta*iv + delta*ivtreat + beta[1]*x1 + beta[2]*x2 + beta[3]*x3 + u
		d <- data.frame(id, y, iv, treat, ivtreat, x1, x2, x3)
		that <- rep(NA, q1)
		for(k in 1:q1) {
			that[k] <- coef(lm(y ~ iv + ivtreat + x1 + x2 + x3, data = d[d$id == k, ] ))[3]
		}
		res[i] <- CRStest(that)
	}
	
	totres[h] <- mean(res)
	write.table(totres, output)
}
